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We calculate the asymptotic behaviour of the one body density matrix of one-dimensional im- 
penetrable bosons in finite size geometries. Our approach is based on a modification of the Replica 
Method from the theory of disordered systems. We obtain explicit expressions for oscillating terms, 
similar to fermionic Friedel oscillations. These terms are universal and originate from the strong 
D ' short-range correlations between bosons in one dimension. 
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I. INTRODUCTION 



Recently, one-dimensional (ID) Bose gases have been created in long cylindrical traps by tightly confining the 
transverse motion of particles to zero-point oscillations pj . These experiments revived an interest in exactly solvable 
one dimensional models of statistical physics, in particular in the Lieb-Liniger model 0] for one-dimensional bosons 
interacting via a delta-function potential. In this case the Bethe Ansatz solution accounts for the ground state 
properties, spectrum of elementary excitations Q and thermodynamics Q| . In contrast to these properties, following 
from the solution of relatively simple integral equations, the correlation functions are not easily obtained from the 
Bethe Ansatz due to an extremely complicated form of the wave functions. In the limiting cases of weak and strong 
interactions closed analytical results can be found as perturbative expansions 0, @ . 
' Another drawback of the Bethe Ansatz is the requirement of periodic boundary conditions and thermodynamic 
limit. This precludes the study of finite size effects, except for a number of cases in which a special symmetry of the 
confining potential Q is available. In most cases 0, E| the finite geometry was treated relying on the local density 
O , approximation, where the macroscopic length scale induced by the confining potential is assumed to be well separated 

from the microscopic correlation length emerging from the Lieb-Liniger solution. 
t^J- | A special case that allows one to go beyond the local density approximation, is the limiting case of an infinite 
£> ■ coupling constant. Then the interactions are taken into account by mapping the system onto free fermions as shown 
by Girardeau [T^j. The resulting system of impenetrable bosons, or "Tonks gas" shares many of its properties 
with free fermions. In the present stage, achieving this strong-coupling Tonks- Girardeau regime is one of the main 
experimental goals [Til IT2I fl3| . In this respect, it is desirable to have results for finite systems beyond the local density 
approximation and to compare them directly with experimental data. 

From a theoretical point of view, the fermionic mapping leads to a considerable simplification of the general Bethe 
Ansatz expression for the wave functions. Historically, even before the Bethe Ansatz solution became available, 
■ Girardeau ^(| was the first to note the boson-fermion correspondence. He found a simple expression for the ground 
state wave function of N bosons in arbitrary potential V(x) in the form of an absolute value of the fermionic Slater 
£h , determinant: 

32! ®(xi,...,xn) = detk,i[<Pk(xi)] , (1) 

P . 

o . 

q , where ipk(x) are one particle wave functions in the potential V(x). 

The equivalence of impenetrable bosons and free fermions can be stated as the equivalence of correlation functions. 
. i-H ! Indeed, any correlation function of the density is given by the corresponding expression for fermions, since it is 
diagonal in field operators and does not involve phase correlations. The same holds for the ground state energy 
and the spectrum of elementary excitations. However, off-diagonal correlation functions of impenetrable bosons are 
different from those of fermions, due to the presence of the absolute value in Q. This drastically changes phase 
correlations. The simplest and well studied example of such off-diagonal correlation function is the one-body density 
matrix 

gi(t,t') = N I dx 2 . . .dx N $*(t,x 2 , . . . ,x N )$(t' ,x 2 , . . . ,x N ), (2) 
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This quantity is of major importance for bosonic systems, since the eigenvalues of g\ show the presence or absence 
of Bose-Einstein condensation according to the criterion of Penrose and Onsager [l4| . In the translationally invariant 
case, gi depends only on the relative distance x = t — t' and its Fourier transform with respect to x is the momentum 
distribution of particles in the ground state. In this case the condensation would manifest itself as a macroscopic 
occupation of the zero momentum state. 

The problem of calculating the one-body density matrix, or equivalently the momentum distribution of impenetrable 
bosons, has a long history in mathematical physics. First considered in 1963 by T.D. Schultz |15| . the one-body density 
matrix was found in the form of a determinant with special properties, that is a Toeplitz determinant |16|. Using the 
known asymptotics of the Toeplitz determinants it was possible to prove the absence of Bose-Einstein condensation 
by showing the power-law decay of the one-body density matrix at large distance. The precise form of this power law 
was obtained later by Lenard in [l^. His calculation resulted in the following long distance behaviour: 

flip) = Poo ,„., 
n \k F x\V 2 ' { ' 

where n — kp/it is the density of particles. In their tour de force, Vaidya and Tracy |l8j| calculated from the first 
principles the asymptotic long-distance behaviour of the one body density matrix (the short distance behaviour has 
also been calculated) and found the expression in as the leading term with p^ = Tre 1 / 2 2~ 1 / 3 A~ (i = 0.92418. . . 
where A — exp (1/12 — £'(— 1)) = 1.2824 ... is Glaisher's constant, related to the Riemann zeta function ((z). They 
also succeeded in obtaining the sub-leading terms. This work has been extended to higher order terms by Jimbo et al. 
|19| who related g\ (x) to the solution of a certain non- linear differential equation, the Painleve equation of the fifth 
kind. The general structure of the large-distance expansion of the one-body density matrix consists of trigonometric 
functions, cosines or sines of multiples of 2kpx, each such trigonometric term is multiplied by series of even or odd 
powers of 1/kpx, respectively. Using a hydrodynamic approach, Haldane poj has shown that this structure is a 
general property of any one-dimensional compressible liquid, bosonic or fermionic. This method is unable to predict 
exact coefficients of the 2kpx harmonics, which should be calculated using exact methods. For example, according to 
[Tsl IT^ | the first oscillatory correction to Eq. J3J is given by 

qi(x) lcos2fci?a: ,,. 

a K ' - 1 = A — (4) 

9i(x) + 8 {k F xY K 1 

Together with p^ , the coefficients of sub-leading oscillatory terms provide the full information on the long distance 
asymptotics of the one body density matrix. It is then desirable to have these coefficients in a simple analytical form 
or to be able to calculate them by a perturbation theory. 

Experimental conditions for obtaining the Tonks- Girardeau regime require a small number of particles, which in an 
isolated system can be as small as N ~ 100 atoms 12]. It is then important to be able to extend the results of Eqs. (0 
and to finite size geometries. The case of harmonic confinement is directly related to experiments. It was recently 
studied analytically in 0| by using Haldane's hydrodynamic approach and, independently, in |2lj | by the Coulomb 
gas analogy (for numerical results, see the work |22| and references therein). The expression for the leading term in 
the one-body density matrix was obtained, generalising the expression in Eq. © to a non-uniform density profile. 
The case of a circular geometry has been considered by Lenard |23j | who conjectured the main smooth contribution, 
analogous to J3J in the form: 

- f \ Poo ( r\ 

9n a > = TW7- R72 ' ( 5 

\N sin ira\ L '' ! 

where 2ira is the angle between two points on the circle. This result was justified rigorously by Widom |24j using 
the theory of Toeplitz determinants. It is worth mentioning here that this result follows straightforwardly from the 
conformal field theory (see e.g. [2{| ), which to some extent is equivalent to the hydrodynamic pfjfand Coulomb gas 0, 
|2(J approaches. These methods are capable to predict long-wavelength, large scale behaviour of correlation functions, 
while giving only qualitative answers for details on a scale of mean inter-particle distance. In both geometries, 
harmonic and circular, the finite size corrections, analogous to the expression in Eq. remain unknown. 

Here we present the first calculation of finite size corrections to the one-body density matrix for both harmonic 
confinement and circular geometry. We compare our analytical expressions with numerical calculations based on the 
exact representation of the one-body density matrix as a Toeplitz determinant and find excellent agreement. We 
also take on the task of reproducing Vaidya and Tracy asymptotic expansion in the thermodynamic limit using both 
systems as a starting point and increasing the number of particle and the system size in a way to preserve a constant 
density. Our calculations are the first to resolve analytically the sign ambiguity of the oscillating terms, similar to 
that in Eq. which appears in the studies of Painleve representation for the one-body density matrix [lJj, |27j . We 
find that the sign of all oscillating terms, such as J3J in the Vaidya and Tracy expansion [lg should be reversed |28| . 
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Our calculations are based on a modification of the Replica Method developed in the theory of disordered systems 
[2^1 and applied recently to random matrices [3(1 0^1 03 and Calogero-Sutherland models |33|. Recent progress 
|34l Isfij in calculations based on the Replica Method has elucidated its intimate relation to the theory of non-linear 
integrable systems. The present work gives yet another example of this interconnection. 

The paper is organised as follows. In Section[n]we present our method and illustrate it by calculating the amplitude 
of removing a particle from the ground state of harmonically trapped impenetrable bosons. The one-body density 
matrix for harmonically confined bosons is calculated in Section IIHI Section IIVI describes the calculation of the one 
body density matrix in the circular geometry. We calculate sub-leading terms using perturbation theory in SectionlVl 
The conclusions are presented in Section IVII and mathematical details are given in two Appendices. 



II. DESCRIPTION OF THE METHOD 



The task of reproducing and extending the original calculations of Vaidya and Tracy to finite systems is obscured 
by the technical complexity of their method which consist in the asymptotic expansion of a Toeplitz determinant 
assuming the size of the determinant tending to infinity in order to reproduce the continuous limit. It is therefore 
important to have an alternative way of representing the one-body density matrix. One starts with the representation 
of the one-body density matrix as an iV-dimensional integral , first proposed by Lenard |17| : 



1 r 1 \ 
gi (a) = ± / d N 6\A N (^\...,^ K )\ ]l\l-^\\^-^'\ = (^]l\l-zi\\^-zi\) (6) 
' ^° {=1 \i=i I 

which follows immediately from the definition J2J. We have chosen here a circular geometry of N + 1 particles 
described by cyclic coordinates < 2tt6i < 2tt with periodic boundary conditions. The ground state wave function 
of N + 1 particles is given by the absolute value of the Slater determinant (JTJ composed of plane waves (pk+i(zf) = 
exp(2irik9i) = z, . It is then identified with the absolute value of the Vandermonde determinant 

deti<k,i<N+i[<Pk(zi)]=AN+t(z) = A(zi,...,zii+i)= Yl ( z i- z i)- ( 7 ) 

l<j<fc<W+l 

The above expression is factorised straightforwardly to yield the expression being integrated in 10. Consider now a 
positive integer n and correlation function 

Z 2n (a) = (]J(l-z l f n (e 2 ™-z l ) 2n Y (8) 

Our method is based on the fact that the one-body density matrix can be obtained from Z2 n by suitable analytical 
continuation to n = 1/2. It happens that Z 2n can be evaluated straightforwardly in the asymptotic large N limit. In 
this respect the expression (JHJ supplemented with proper procedure for analytic continuation n —> 1/2 is the desired 
alternative representation of the one-body density matrix. 

The idea of calculating the averages of the absolute value of a non-positive definite function was put forward by 
Kurchan |36| and was named a modification of the Replica trick. In the present context replica means the following. 
In [23, l3ll l32l l33l ] the correlation function (JBJ) was expressed using a dual representation involving an integral over 
components of a n-dimensional field, with the action symmetric under rotations in the space of the components. This 
zero-dimensional field theory was considered in the n — ► limit to obtain density-density correlation function. In the 
present work a different limit n — > 1/2 is taken to reproduce the off-diagonal correlation function J()J). However the 
idea of the analytical continuation in n is common to the above mentioned works and we use heavily the techniques 
introduced in |3l| . For example, we use the same dual representation of JBJ and evaluate it in the asymptotic large 
N limit. 

For n integer, the result consists of a main smooth contribution and exactly n oscillatory corrections. This is an 
expected structure of one-dimensional correlation functions conjectured by Haldane |2f)j | form his hydrodynamical 
approach. The sensible analytic continuation n — ► 1/2 is then performed in a way to preserve this structure and the 
resulting expression is believed to represent the large N limit of the (unknown) analytical continuation in n. The lack 
of rigour in this approach is shared by most replica calculations and is justified a posteriori by remarkably transparent 
resulting expressions which are in full agreement with the numerics as well as the known analytical results. 

To demonstrate the method in detail and set up notations we first calculate the ground state amplitude A(t) of 
impenetrable bosons in a harmonic potential: 

A{t) = N {V(t)) N+1 (9) 
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This quantity describes the probability amplitude to remove a particle at position t by acting with annihilation 
operator ^(t) from the ground state of N + 1 particles and leave the system in the ground state of N particles. It can 
be considered as a many-body wave function of the removed particle. We consider a geometry different from circular, 
since for the latter the ground state amplitude is just square root of mean density independently of the position due to 
the translational invariance. We consider the system confined by harmonic potential V(x) — mu> 2 x 2 /2 for which A(t) 
has a non-trivial position dependence. The one particle orbitals are given by eigenfunctions of harmonic oscillator: 

",„(.<■) = -=U ni \ .r\/- ) , --<-. ( ^) 2 2 m ml, (10) 



and H m are Hermite polynomials. We measure the distances in units of half of the Thomas-Fermi radius R — 
2hN /muj, corresponding to half the size of the particles cloud in the large N limit. The ground state wave function 
is obtained similarly to the uniform case taking the absolute value of the fermionic Slater determinant: 




$(x\,...,x N ) = det ; . m [0 m _i(x ; )] = = A N (x) 



e-T £, i. (11) 



The last identity follows from the linearity of the determinant with respect to its columns, which enables us to write 
the determinant of Hermite polynomials H m (x) as a determinant of their leading monomials H m (x) ~ x m . The 
normalisation constant is given by Selberg integral |37l l38j of Hermite type 

S N (X)= AA^( 3 ;)e-^^A- w2 / 2 (2i) A '/ 2 I]r(l+j). (12) 

J-oo j=1 

Using the definition JjJ) of the ground state amplitude leads to the expression 

A(f) - vmrnm L A a ~ fe)e ~ f E * JJ > - * 1 (13) 

To deal with absolute value in this expression we use the identity 

\t- yj \^ = (t- yj ) 2n , (14) 

valid for integer n. Performing the analytical continuation n — > 1/2 in the end of the calculations, one recovers the 
expression (|13|) ■ To treat the ground state amplitude and one-body density matrix on the same footing it is convenient 
to consider a more general quantity 

i „oo N m 

Z m (t u ...,t m ) = —— d^A^e-^^HHiU-yj) (15) 

The variables t ai different in general, arc later taken equal to a single value t, 

Z(t)= lim Z 2n (t)= lim Z 2n (t,...,t) (16) 

n— >l/2 n— >l/2 ^ ^ ^ 

2n 

to recover Eq. Ijl3(l up to a normalisation: 



It is crucial that one cannot set 2n = 1 directly in (|15|) . which means that Z(t) in (|16(l is different from Z\(t). The 
latter is nothing but a fermionic ground state amplitude obtained by removing the absolute value in H13J1 • It is given 
by the wave function 1)10(1 of TV + 1-th particle removed from the system and results in an expression that oscillates 
rapidly around zero with period equal to the mean inter-particle separation. In contrast, the bosonic ground state 
amplitude, obtained in the large N limit, is expected to contain a smooth positive leading term which is a direct 
consequence of positivity of the ground state wave function l|ll|l . Our claim is that it can be obtained from Z(t). 

Recently the leading smooth contribution to Z2 n (t) for integer n was evaluated in the large N limit by Brezin and 
Hikami |39| and it was shown to survive the analytical continuation in n off the integers. This suggests that we deal 



N +2 
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with two functions of variable n, one for fermions and one for bosons, or equivalently one function with two branches, 
a fermionic one and a bosonic one which coincide at integer n but become different as n is moved away from integers. 
The goal of the present calculation is the generalisation of the bosonic analytic continuation considered in |3j| to the 
whole asymptotic expression for the ground state amplitude. 

To illustrate the ideas above we start with a remarkable duality transformation, which represents the TV-fold integral 
Z2n in Eqs. I|15jl. <|16|l as an integral over m — 2n variables: 

1 r°° 

3*(*)=o-?An/ ^xAt n (x)e-»Z a s M (18) 

<J2n(.JV) J_ OQ 

with an effective action 

S(x,t) = ^^-lnx+^. (19) 

This representation is exact and its proof can be found in the mathematical literature (see [4p| | and references therein). 
For the case of harmonic confinement the proof was presented in [3l| using the Random Matrix Theory. This proof 
has the advantage to be readily extended to deal with two-point correlation functions in harmonic potential. In 
Appendix ^ we present yet another proof using second quantization for fermions. 

The dual representation (|18|l provides us with an alternative representation of the TV-fold integral (|13[1 and is the 
starting point of the large TV asymptotic expansion. The dominant contribution comes from the saddle points of the 
action fP9")l . situated at 



*± = f ± ] /l-j = ±e ±i *, sin^=i (20) 

for t inside the ground state density support — 2 < t < 2. The stationary value of the action and its second derivative 
at these points read 

e^ 2 ^ ni 

S(x ± ,t) = S ± = —Ti4>±Y> ( 21 ) 
S"(x±,t) = S'± = 2e =F< *cos0 (22) 

To take into account all possible saddle points we put I variables x a in the vicinity of X- and 2n — I variables in the 
vicinity of such that 



x a = X-+£ a /vN, a = l,...,l 

x b = X+ + &/VN, b = l + l,...,2n. (23) 

In what follows the saddle points with I = or 2n will be referred to as replica symmetric for obvious reasons. The 
other saddle points break the replica symmetry to some degree, the maximal degree of replica symmetry breaking 
happens for I = n. We shall see that the analytical continuation of the contribution of this saddle point leads to the 
result of [3|j. On the other hand, in the calculation of the fermionic ground state amplitude, the replica symmetry is 
preserved, since in this case the contribution comes from the replica symmetric saddle points. 

The pre-exponential factor in l|18|) given by Vandermonde determinant vanishes at each saddle point (|23() and should 
be expanded as 

/ , i(i-l)+(2n-2)(2n-i-l) 

AUx) = l^=j (x_ - x + f^Af^ a )Aj n _^ b ) (24) 

to yield a non-zero result. Using this expansion the integration of the fluctuations near each saddle point is performed 
by using the Selberg integral 1121) . Combining the results of integration with the main saddle point contribution and 
summing over all saddle points we represent the original integral (|18|) as a sum of 2n + 1 terms 



Bn. , (r,-r N |2/(2n-i) 

sT) P (V^T) (2 ""° 2 



Z2n(t) = E FlN l ^-D . g± p ^_ {2n _ i)2 e -NlS-N [2 n-l)S + (25) 

The factors 



1=0 



F i (2n\ nU r(a + 1) Eg r(b + i) _ -rr r( a ) 



l J n^rfc+l) l\T(2n+l-a) 
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combine the numerical factors originating from Selberg integral and the binomial coefficient equal to the number of 
ways to choose I variables x a in the vicinity of x_ out of total 2n variables. 

As 1 5+ 1 = 15-1, the N dependence of each term in (|25(l enters only through the factor N l ^ 2n ~ l \ Therefore the 
central term of the sum (|25|) with / = n corresponding to the maximal degree of the replica symmetry breaking gives 
the dominant non-oscillating contribution to Zi n (t). The other terms, including the replica symmetric terms I = 0. 2n 
are at least 1/N times smaller and oscillate rapidly. This is to be contrasted with the replica approach in [sH Issl] 
where in the limit n — > the dominant contribution comes from the replica symmetric points and the role of the 
saddle points with broken replica symmetry is to provide oscillatory corrections. 

Observing the behaviour of l(2"B|) for integer n we are lead to the conjecture that the right analytical continuation off 
the integer n following the bosonic branch would preserve this form: one smooth term plus oscillatory corrections. We 
now prepare the expression l|25|) for the analytical continuation. To this end it is convenient to change the summation 

variable I — k + n and factorise the amplitude (j26(l as F\ n = AnD^, where the first factor is given by the product 

a « - n Tv/i*? v w 

^-^ 1 (2n + 1 — a) 

where n enters explicitly as the range of product. The analytical continuation of this factor to non- integer values of 
n is described, for example in 0,123. We reproduce it in Appendix iBl where it is shown that p^ = A\^ 2 /\/2. 
The analytic continuation of the second factor 

(n)_-rT T(n + a) _ T(n + k) (n) (n) _ („) („) , . 

U k -ll r(n + 1 _ o) - r(n + 1 _ jfe) ^fc-i. fc >°i D o -1, » k ~ D -k l^j 

to non- integer values of n is straightforward. For integer n the coefficients fl" vanish for |fc| > n due to the 
divergence of the gamma function in the denominator, so for integer n the sum in (|25|) can be formally extended to 

(n) 

—oo < k < +oo. For general values of n, the factor D k has a non-zero value for any k, which results in genuine 
infinite series. These series are asymptotic, rather than convergent, but it can be shown by using the Stirling formula 
that the coefficients decrease very fast for small values of k. Hence, for large enough N, few terms around k = 
provide an excellent approximation for the sum. We rewrite the series (12511 using the new summation variable: 



,_ x i.vYros-V,) 



Z 2n {t) = A n {2Ncos(t>) n e - Nncos2 ' 1 ' V * e ^NkQ-2 in k4>^ ( 29 ) 



where O(i) = 2<ft + sin 2(f) + ir — 2ir J^ 2 p(s)ds is expressed as an integral of the mean density of particles p(t) 
(I /it) cos</>, given by the celebrated Wigner semi-circle law: 



p(t) = ±Jl- t l (30) 



The result has the expected form consisting of sum of smooth and oscillating parts: 
Z 2n (t) = A n (27rNf e -^(i-* 2 /2) 



x [p(t)Y 



D (n) 



2 zZ " J? cos (^©(^ + 2nk(f>{t)) 



k=l 



(31) 



This is the desired asymptotic representation of Z n . For an integer n it has a simple general structure: a smooth main 
contribution plus oscillatory corrections. This structure is preserved if analytical continuation n — > 1/2 is applied 
term by term. The validity of the change of summation variable I = n + k for a non- integer n is explained as follows: 
we assume that in fact the infinite asymptotic series (|29(l is the correct expression of the large N asymptotics for the 
analytical continuation of Z 2n (t) to arbitrary n. For integer n the sum terminates and can be rewritten as l)25|) after 
the corresponding change of the summation index. 

Multiplying l|31|) by the normalisation factors defined in l|36|l and using the asymptotic expansion 

S N (N) 1 N N +i 



S N+ i(N) V2^(N + l)l 2ttN ' 
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we obtain the desired asymptotic expression for the ground state amplitude 



A(t) 



D 



(1/2) 



1/2, D { 2 1/2) -- 



1/4 



1 + 2E 



D 



(1/2) 



k=l 



■cos [kNB(t) + kcj>{t)] 



-3/8,. 



D 



(1/2) _ g (l/2) r (fe+|) 



fe+1 



r( 



(33) 
(34) 



The explicit values of the coefficients of oscillatory terms are obtained from (J2SJ. 

In recent papers |34l l35l | the functions Z m (t) were shown to satisfy a remarkable recursion relation, a Toda Lattice 
hierarchy, extensively studied in the theory of non-linear integrable systems |4l| . The same Toda Lattice equations 
relates the solutions of Painleve equations, where m enters as a parameter not restricted to integer values |42j. This 
observation was crucial for exact evaluation of the replica limit m — > without relying on the large N asymptotics. 

In the present framework the quantity of interest is Z m (t) for m — > 1, so according to [3~H l42j | it can be related 
straightforwardly to the solution of Painleve IV equation. The boundary conditions, obtained from the small t 
expansion of the ground state amplitude (|13fl distinguishes between bosonic and fermionic branches of Z\(t). Our 
analytic continuation with broken replica symmetry chooses automatically the correct boundary conditions for bosons. 

To justify numerically the result (|34|l we have calculated the ground state amplitude as a determinant 



A(t) 



N \ e ~Nt 2 /4 



y/S N (N)S N+1 (N) l<j,k<N 



det A j+k - 2 (t), 



(35) 



where the matrix elements can be expressed by gamma function and incomplete gamma function j(a,x) as 



m+l Nt 2 



The results are presented in Fig. for N = 20 particles and — 1 < t < 1. In the central part of the cloud the 
ground state amplitude is indeed well approximated by few terms in the expansion l|34l) . Further improvement can 
be achieved using the perturbation theory around each saddle point similarly to the perturbation theory described in 
Section in ther case of circular geometry. Close to the edges of the atomic cloud, which in our units correspond to 
t = ±2, the saddle point approximation we used to derive l|34(l breaks down since for t approaching one of the edges 
the saddle points l|2(J|l coalesce at ±i and the second derivative of the action l|22|l vanishes. In this case the analysis 
of the higher order expansion of action (1221 is required and is beyond the scope of this work. 

Thus we see that our method of calculation provides the exact large N asymptotic of the ground state amplitude. 
We extend it to the calculation of the one-body density matrix in the subsequent sections for the case of particles in 
harmonic potential and circular geometry. 



III. ONE-BODY DENSITY MATRIX IN HARMONIC POTENTIAL 



The one-body density matrix can be written using the definition @ together with the expression JQ| for the ground 
state function 

9 1 (t,t') = (N + l)^Le-^^Z(t,t'), (36) 

where Z(t,t') is obtained from (|15|l according to the following rule: 

Z(t,t')= lim Z 4n (t,t')= lim Z in (t, . . .t,t' , . . A 1 ) (37) 

2n 2n 

It is again important to calculate Z± n (t,t') in the large N limit before taking analytical continuation n — ► 1/2, 
otherwise the result will be just the one-body density matrix for one-dimensional fermions. 
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FIG. 1: Comparison between asymptotic expression for the ground state amplitude 13411 (solid line) of N — 20 particles in 
harmonic potential evaluated up to the k = 1 terms and exact results (dotted line) based on numerical evaluation of the 
determinant 13511 . The contribution given by the smooth term k = in 134H is represented by dashed line. The ground state 
amplitude is normalised in a way that the smooth part equals unity for t — 0. 



The duality transformation can be worked out in the case of two variables as explained in the Appendix ^ and 
yields 

^n(M') = ^AFT r rf2 "£ d2Hx> A 2n ( x ) A tn ( X> ) ^ a ' = 1 ^" \ ~ N T. S(x a ,t) £ -N Y, S«, ,t') ^ (3g) 

°2nK ly ) J-oo [l (t — t')\ 

with the same definition l|19(l for the effective action, S(x,t). 

Before calculating the one-body density matrix by the saddle point method we would like to remark that this 
correlation function depends on the scaling of the distance t — t' when TV goes to infinity. We distinguish two limits: 
macroscopic, when t,t' are of order of the cloud size and microscopic, or, more precisely, mesoscopic limit, when t, t' 
remains finite on the scale of mean inter-particle distance 1/N pit) equal to ir/N in the centre of the cloud. The 
second limit is called mesoscopic, since though t, t' are small on the scale of the cloud size, we are interested in the 
asymptotic behaviour of correlations t — t'^> 1/N, so there are still a very large number of particles between t and t 1 . 

In the macroscopic limit, the variables are changed in the following way: 

x a = x-+£ a /VN, a = l,...,n + k 

xt = x+ + vJV, b = n+ k + 1, . . . ,2n 

<• = x' + +? a ,/VN, a' = l,...,n + k' 

x' v = x'_+£,' b ,/VN, b' = n + k' + l,...,2n (39) 

where x± = x±it) and x± — x±it') defined in Eq. 1)20(1 . The stationary value of the action and fluctuation integrals 
are done exactly as in the last section. The only new factor is the double product in the integral, calculated at the 
saddle point (k,k') with the result 

kk' 

e ~2m(fc0-fc'0') MO) 



This factor represents the interaction between saddle points in x a and x' b and as we shall see it is crucial for obtaining 
the correct expression for g\ in the mesoscopic limit. Combining this factor with each (fc, k') saddle point contribution 



2n In 



n n(^-<')=* 4n2 i 2n(fe+fe,) i*-^T 



a 4>+<t> 



a=l a' = l 



Sill 
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and multiplying by the number of ways to distribute x, x' among different stationary values we get the integral l|38|) 
as a double sum 



V (4 +>\ 42 /o AT\ 2n -Nn(2- t 2 +t' 2 )/2) „ [p(t)p(t ) 

Z in (t,t ) = A n (2ttA) e V V )i ) x — 

f - t'\ 



2n 2 



E E 

fc— — 00 k'=—oo 



<t>+¥ 



sin ■ 



2U'' 



• -Tfe 2 
8Nir 3 p 3 (t) 



-4Nk&— Hnk$ _ 



iNk &' +Aink' </>' 



(41) 



where the summations are extended to infinity, relying on the factors which cut off the finite number of terms 
in the sum. 

We see immediately that again the most replica asymmetric saddle point (k, k') = in each set of variables x a and 
x' a , provides the dominant smooth contribution, which after analytic continuation n — > 1/2 and normalisation given 
in l|3tj|) yields the macroscopic one-body density matrix 



9i (*,*') = Poo — 



NY WW)]' 



t-t' 



(42) 



This result is identical to that of Ref . [2lJ and agrees completely with the functional form deduced in . The finite 
size correction to <|42[) are obtained by taking n = 1/2 in each terms summed up in (|41|l with the result 



9l(t,f) 1 _ V-> r _ 1 x(fc+fe')/2 n (l/2) n (l/2) 

9i(t,t') ^ [ ' k k ' 

(fe,fc')#(0,0) 



4>+<P' 



2kk' 



e -iNk&-2ik4>+iNk'0'+2ik'4>' 

n k 2 r , k' 2 

8Nn 3 p 3 (t) 8NTr 3 p 3 (t) 



(43) 



Going to the mesoscopic limit we focus on the centre of the potential t + 1' = and define the scaling variable x 
such that t — t' = x/N. The factor in (|4*T)l becomes 



cos 



4>+<P' 



sin ■ 



2kk' 



AN 

X 



2kk' 



(44) 



Adding powers of TV in this expression with those in (|41|l we see that the diagonal elements k = k' in the sum are of 
order 1, while the off-diagonal elements are at least smaller by factor 1/N. Therefore we anticipate that the mesoscopic 
limit is given by the diagonal part of the sum (|41() . 

In principle, one should reconsider the asymptotic expansion l|41(l in the mesoscopic limit, since in this case the 
distance between saddle point is of order 1/N, which can change the contribution of the fluctuations. These calculations 
can be done straightforwardly along the lines of Ref. ■ Somewhat surprisingly, the conclusion is that the asymptotic 
limit N — > 00, and the mesoscopic limit t — i' — > commute, so one can change the variables to x directly in the sum 
gT}. Apart from the factor |@IJ| the other factors simplify as p(t) = p{t') — ► 1/tt, 6 — 9' = 4<p — 4<fi' = 2x/N and one 
gets 



Z An (x) = Ale- N < 2 -( t2+t ' 2 )' 2 



4n A r 



! + 2 E(-i 



2nk 



fc=l 



D 



in) 



2 cos 2kx 
(2x) 2k2 



Normalising and putting n = 1/2 one recovers the dominant term: 

Poo 



l N 

01 W = — 



7T / X 



1/2 



(45) 



(46) 



which is identical to the expression © of Vaidya and Tracy if one identifies kp — nN p(0) = N and normalises to the 
density Np(0) = N/tt in the centre of the cloud. By the same procedure we obtain the oscillatory corrections 



9i( x ) 



1 = 2E(-!) 



fe=i 



D 



(1/2) 



2 cos2fca; 



,2k 2 



(2x 



(47) 
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Using the explicit value D^ 2 ^ = 1/2 we see that the k = 1 term is identical to the the leading term (0J of Vaidya 
and Tracy expansion [18j apart from the sign. To our knowledge it is the first analytical justification of the sign 
inconsistency in the Vaidya and Tracy asymptotic expression first noticed in \27\ by using numerical solution of 
Painleve equation. The same sign change appears in the thermodynamical limit of circular geometry which we 
consider in the next section. 

The expression l|43|) and Ij47(l together with the expressions l|34|l of the coefficients D^ 2 ^ provide only the leading 
non-perturbative contribution of each saddle point and do not contain the sub-leading terms, which arise, for example, 
from the deviation of the action from its second order Taylor expansion near saddle point or additional terms in the 
large TV expansion of the double product H44J1 • These terms can be treated by a perturbation theory near each saddle 
point and we show in Section how such perturbation theory can be constructed for the case of circular geometry. 



IV. ONE-BODY DENSITY MATRIX IN CIRCULAR GEOMETRY 



Our method is also able to give the one-body density matrix in the circular geometry, i. e. for N + 1 particles on a 
ring of length L with periodic boundary conditions. The one-body density matrix is given by the expression @. The 
corrections to the leading smooth term ® in the long-distance expansion of the one-body density matrix remained 
unknown to the best of our knowledge (see discussion in [2l]]). We now show how they can be obtained with our 
method. Let us consider the quantity (JHJ) rewritten explicitly as an iV-dimensional integral 

Z 2n (t) = - 1 / 1/2 d N 8 | A* (e 2 - e )| 2 n 1 1 + ^ lBl T I* + ^ ^ > (48) 

M N (2n,2n, 1) J_ 1/2 fj^ 

where t — exp(27ria) and we have normalised Zm(X) — 1j introducing the normalisation constant Mjv(2n, 2n, 1). Its 
value is given by Morris integral of Random Matrix Theory [3^, but its precise value is not needed. Taking n — ► 1/2 
in Zm gives the density matrix normalised to the density gi/n. Now we write 

|1 + e 2 ^ 9 '\ 2n = e - 2 ™ ie > (1 + e 2 ^ 9 ') 2 " (49) 
and change the integration variables 9i — * + a to obtain 

j.— Nn i-l/2 N 

Z2n{t) ~ M N (2n,2n,l)J_ 1/2 d )Hl e ' 1 + 6 ' ^ + > " (5 ° } 



One observes that the transformation from (|48() to the last expression is possible only when n is integer, otherwise the 
change of variables is not permitted due to discontinuity of the phase in l|49|l . It is parallel to the representation (|L4|) 
we used for the harmonic confinement. As in the last section we proceed with the representation (| 5 Of) and assume 
that it remains valid for any n. 

The integral H50JI has a remarkable dual representation (see eq. (3.41) in 40] ) by an integral over n variables 



f -Nn i-l 2 ™ 

Z2n(t)= „ , nnn / <F n xAt n (x)l[(l-(l-t)x a ) N , (51) 
J2n(U, U, I) Jo 



where the normalisation constant is given by Selberg integral 



5^(0,0,1) = / d 2n xl±l n {x) = J] (52) 
Jo „ = i r ( 2 « + a ) 

In the right hand side of l|51|l the number of particles TV appears only as a parameter. This representation is a direct 
analogy of i|18|l and (|38fl and allows us to obtain the asymptotic expression for Z 2n suitable for analytic continuation 
in n. In the large N limit the integrand in (|51|l oscillates rapidly and the main contribution comes from the endpoints 
x± = 1,0 which are the only stationary points of the phase. We change variables near each endpoint 

a — !,...,{ 



N(\-ty 



x b = x + ~ N ( l _ t -iy b = l + l,...,2n, (53) 
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The integrand in i|51|) simplifies in the large N limit: 



(l-(l-i)O 



N 



e a = 1, . . . , I 

t N e-t°, a = l + l,.. 



,2n 



and the integration measure including the Vandermonde determinant is factorised as 

l 2 / -, x (2n-2) 2 



1 



1 



^N(l-t)J \N(l-t- 1 ) / 
The remaining integrals are calculated using Laguerre variant of the Selberg formula 



f°° . m . J 

J,(A) = / dftaAftta) J] e- A? " = A-' 2 [] r(a)r(l + a) 

J ° o=l a=l 



(54) 



(55) 



(56) 



Multiplying the contribution of each saddle point by the number of ways to distribute variables we obtain the asymp- 
totics of the integral (|51() as a sum of 2n + 1 terms: 



2 ™ [rpl ] 2 f (N+2n)(n-l) 

z 2n (t) = n^ + e) E(-D 2 " ( "" ) 1 r,"U + ( 2 n. 



(2X)* 2 +(2«-0 2 



(57) 



where we have introduced X = N sin 7ra and the factors F\ n are defined in l|26|) . We note that due to the translation 
invariance the resulting expansion is given by simple sum over saddle points and not a double sum as in the case l|41l) 
of harmonic potential. 

Changing the summation variable to k = I — n and factorising the amplitudes F\ n according to the definitions l|27Jl , 
(|28|l we obtain finally 



(27V | sin7ra|)^ n \ 



2 cos [2kn(N + 2n)a] 
(4/V 2 sin 2 7ra) fe2 



(58) 



Now we are in a position to take the limit n — > 1/2 which results in the desired corrections to the smooth part JSJ) of 
the one-body density matrix: 



gi(oQ 



-l = 2^(-l) fc 



(1/2) 



k=l 



2 cos2fc7r(iV + l)a 
(4/V 2 sin 2 ?ra) 



(59) 



In the thermodynamic limit Nira — Rfx, when N — ► oo we recover the result (I47|l . Again we note the sign difference 
between the k = 1 term of l|59|l and that of Vaidya and Tracy l@J. As in the case of harmonic confinement only 
the leading contribution of each saddle point is included in the expansion l|59|l . The circular geometry is particularly 
suitable for discussing the perturbative corrections, which are calculated in the next section. 



V. PERTURBATION THEORY 



Up to now we have considered only the leading non-perturbative contribution of each saddle point in the expansion 
(|58|l of Z2n- To provide all the corrections to a given order in 1/N we have to deal with the sub-leading terms by 
perturbation theory. This amounts to multiplying the contribution of each saddle point in l|57|l by the following 
average 



of the function 



n-\-k n — k 




b=l 



(60) 



/ n-k / f ,2 >,3 \ \ 
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where we have defined £' b = 6+n+fc for b = 1, . . . , n — k and 



A 



(ab) 



I 

2 ,f 



Triot —ma 
— e 



6 



(62) 



In l|61l) the product comes from the neglected terms in the factorisation of the Vandermonde determinant, while the 
exponents represent the corrections to the leading term l|54|) in the 1 /N expansion of the action. 

The perturbation theory consists in expanding F(£, £') up to the desired order in 1/N (recall that 1/X = l/7Vsin7ra) 
and averaging (|(jl)f) term by term with the unperturbed action of each saddle point k. It follows from the very essence 
of our method, reflected in the structure of (|SU|) . that the sets of variables 6 and f£ are independent, so that their 
averages factorise 



/(o$(0) = (/(o) (<?(£'; 



i—k 



\ Ira l m (L) Jo ^ 



(63) 



where the subscript reflects the number of integration variables in each factor. The remaining averages are performed 
using known results from the theory of Selberg integrals (see chapter 17 of the Ref. j^]}). In the following we shall 
need the following results: 



6 



66 



2m 2 , 



m(m — 1), 



666 



to(5to 2 + 1), 
m(m — l)(m — 2), 



Ufa) =m(m-l)(2m-l), 



6 2 66 



6666 



to (8m + 15to 2 - 2m + 3), 
2m(m- 1) 2 (to- 2), 
to(to — 1)(to — 2)(to — 3) 



6 3 6 



66 



to(to — l)(5m 2 — 4m + 3), 
■■ m(m — 1)(2to — l) 2 , 



(64) 



In the first order we have 
2 

X 



n-\-k n — k 



a=l 6=1 



i —k 



6 



n+A; 



n — k 



_^ it-rn, ^ n rv 

_ 2iV (^ a )„_ Li . _ 9aF X! 



+fc 2N 



b=l 



n — k 



— (n + k){n - k) \(n - k)e ma - (n + k)t 



-[(n + k) 3 + (n-k) 3 } 



(65) 



where we have used the fact that the averages are independent of the index a or b and the corresponding sums yield 
a factor n + k or n — k respectively. The calculation of the next order correction to this term is performed similarly 
to the first order. One only has to pay attention to the appearance of identical indexes in the sums. For instance we 
have a sum of diagonal terms: 



X 2 



ab 



(n + k){n — k) 
4X2 



I n-k \ I n+k \ I n+k 



I i 



— k 



(66) 



and non-diagonal ones: 



— ]T (A Qb A a 



v = 



ab^a'b' 



(n + k)(n — k) 
2X~ 2 



(n + k-l)U[ 



+ (n + k)(n-k-l){&& 

-k \ I n — k 



13 



(n-fc-l)(g) +(n-k){n + k-l)U 1 b) 



+ 2((n + fc)(n-*)-l)(&) (#) 

\ / n+fc \ / 1 



-As 



(67) 



Restricting ourselves to the terms up to the second order in perturbation theory we see that only the case k = 
has to be considered, since the leading contribution of k ^ saddle points is already at least of second order in 1/X. 
In this case the first order terms (|65|l simplify to — 4n 3 /A, while calculating the contribution of all the second order 
terms in (|61H . including that of l|66|) and l|67[l yields 



me. 



2X 2 



N 2 



8n q 



23 



-rf - 2n 



(68) 



In every order in perturbation theory the coefficient is given by a polynomial in n, so its analytical continuation to 
n = 1/2 is straightforward. A peculiar feature, familiar in the replica method, of such perturbation theory of this 
n = 1/2 component field is that averages of positive quantities £, £' vanish or become negative. 

Adding up the second order contribution 1|68|) with that of first order 116511 and setting n = 1/2 we obtain the 
factor multiplying the contribution of the k = saddle point in the expansion (JSSJ. Combining it with the leading 
contribution of the k = ±1 saddle point we get the finite size correction up to the second order: 



9i (a) = 



|iVsin na] 1 / 2 



1 



1 

2N 



13 



1 



cos2?r(A^ + l)a 
32N 2 ~ 32N 2 sin 2 na 8N 2 s in 2 ira 



(69) 



We have compared this result with the exact calculation based on a numerical evaluation of the Toeplitz determinant 
representation [TtI l2l| of 51(a). The result of the comparison are presented in Fig. for N = 20 particles. The 
agreement of the two expressions is remarkable. 



-0.02 



3^ 
105 

5 




a 

FIG. 2: Comparison between asymptotic expression 1)69^ for the one-body density matrix (solid line) of iV = 20 particles in 
circular geometry and exact results (dotted line) based on numerical evaluation of Toeplitz determinant. The correction due 
to the finite number of particles C = — 1/2JV + 13/32JV 2 = —0.023984375 is shown and the dashed line represents the mean 
contribution of the oscillatory terms. 

Higher order corrections can be calculated in similar way. However the number of averages proliferates quickly 
and the computations become very tedious. The simplification occurs in the thermodynamic limit, where terms 
proportional to just inverse powers of in (|61|) can be neglected. In this limit only the double product contributes 
in (|61|l and the exponential factors in (|62|l can be set to one from the very beginning. In this case we were able to 
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proceed up to the terms of order 1/X 4 , reproducing the results of the asymptotic expansion of Vaidya and Tracy. 
Identifying in the thermodynamic limit X = N sunra = Nira = k F x we get 



gi{k F x) 



\k F x\V 2 



1 



1 cos(2fci?a:) 3 sm(2k F x) 



33 



1 



93 cos(2k F x) 



1 1 

32 (k F x) 2 8 (k F x) 2 16 (k F x) 3 ' 2048 {k F x) A ' 256 (k F x) 4 



(70) 



To compare this result with the numerics we use the representation for g\ given in the work of Schultz |l5j as 
a continuum limit of the one-body density matrix on a lattice which was a starting point for Vaidya and Tracy 
calculations |l8j| . It allows for the direct study of <?i in the thermodynamical limit if the distance on the lattice d is 
related to the dimensionless distance in l|70l) as nisd — k F x. The filling factor v approaches zero while k F x is held 
fixed in order to reproduce the continuous limit. The results of the comparison are shown in Fig. [3] and the agreement 
is good order by order. In particular, our result (|TC)|> implies that the sign of all the trigonometric terms in Vaidya 
and Tracy expansion should be reversed as in <|T()|> . while the values of numerical coefficients remain unaffected. 




k F x 




k F x 



FIG. 3: Comparison between different orders in l/k F x in the asymptotic expansion I70H and results based on numerical calcula- 
tion of the one-body density matrix g\{k F x) on the lattice with filling factor v — 0.005. a) Dotted line: the Toeplitz determinant 
evaluation of gi(k F x), solid line: gi(k F x) = Poo/j^F^I 1,72 . b) Dotted line: C?' 1 ' (k F x) = (k F x) 2 (gi(k F x)/gi(k F x) — 1), solid line: 
G<$ = -l/32-cos2fc F x/8. c) Dotted line: G (2) (k F x) = k F x{G {1) (k F x) -G$ (k F x)), solid line: G^ (k F x) = —(3/16) sin 2k F x. 
d) Dotted line: G (3) (k F x) = k F x(G (2) (k F x) - G^(k F x)), solid line: G (3) (k F x) = 33/2048 + (93/256) cos2fc F a;. 



VI. CONCLUSIONS 



The main results of our calculations are expressions for the one-body density matrix accounting for finite size effects 
in the cases of harmonic confinement and circular geometry. These expressions consist of a sum of the leading smooth 
term and sub- leading oscillating terms. These terms involve universal amplitudes D^ 2 ^ given explicitly by Eq. 13411 . 
The oscillations occur with periods 2k F , ik F , . . .and are similar to Friedel oscillations in fermionic systems. Their 
physical origin lies in particle correlations on scales of mean inter-particle separation and they represent a generic 
property of all strongly correlated one-dimensional systems. For fermions their existence is due to the Pauli principle 
and is not restricted to ID, while for bosons the Friedel oscillations are due to the effective Pauli principle induced by 
strong interactions, a situation that is possible only in the one-dimensional world. Our results are in agreement with 
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the general structure of the asymptotic expansion of g\ given in [18| and with Haldane's hydrodynamic approach [20j 
or conformal field theory [25| . 

Using semiclassical methods, such as local density approximation, this universal behaviour of correlations can be 
extended for a wide class of external confining potentials with sufficiently smooth behaviour. Our calculations provide 
ab initio justification of the local density approximation in the bulk of bosonic cloud in a harmonic trap which is 
important to current experiments with cold atoms. For this system our method is capable to consider the edges of 
thermodynamical density profile, thus going beyond the local density approximation. Another system of experimental 
interest is the lattice model of impenetrable bosons, produced recently [bl by confining atoms in optical lattices. Given 
an exact expression for the one-body density matrix in terms of Toeplitz determinants |l5l Il7| it is highly desirable 
to obtain its asymptotics without going to the continuum limit, thus revealing the interplay between the short-range 
correlations and the effects of lattice. 

The results were obtained using a novel modification of the replica method which in the context of exactly solvable 
models consists of alternative representation of correlation functions. The phenomenon of replica symmetry breaking 
serves here as a tool to single out the bosonic branch of the one-body density matrix out of various possible analytic 
continuations in the replica index n. It can be considered as a fresh insight in quantum statistics in one dimensions |43j . 
a question to explore in the future. The obvious candidate for this study is the Calogero-Sutherland model. Indeed, 
the ground state wave-function of the Calogero-Sutherland model is proportional to the Vandermonde determinant 
taken to a power A, which characterises the statistical interactions between the particles 0]. Based on equivalence 
of Calogero-Sutherland models and Random Matrices for particular values of statistical interaction parameter A, the 
replica method also provides an efficient way to study correlations of spectral determinants of random matrices directly 
related to averages of the form (jSJ- These objects play an important role in statistical physics, mathematical physics 
and modern combinatorics fsee |39l l45j and references therein). A more distant, but certainly tempting perspective 
is the application of the Replica Method to other integrable models and calculation of their correlation properties. 
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APPENDIX A: DERIVATION OF THE DUALITY FORMULA 



Using the obvious identity among the Vandermonde determinants 

N m 



nTT/, «. ^■N+ m {tx,-..,t m ,y 1 ,...,y N ) 



the expression l|15|l is represented in the second quantisation 



Z m (h, t m ) - t/ ^±!!i e f £ tl _!_ N (^{t 1 )^{t 2 ) . . . ^{t m )) N+m , (A2) 

using the fermionic annihilation operators ip{x) acting between ground state of N and N + m fermions. One uses 
then Wick theorem to calculate this matrix element: 

jv(^(*i)V'(*2) ■ ■•V'(*m))iV+m = det [^jy+fe-i (ij)] • (A3) 

k,l 

The resulting determinant in the right hand side of l|A3(l is constructed using the one-particle wave functions l|l(J|l . 
which have an integral representation 
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Representing in this way the one-particle wave functions in the expansion of the determinant (|A3|) we have the 
following result: 



H) 



Nm 



d m xe~^^ 



(x a -it a ) 2 A-m(xii ■ ■ ■ ; x m) 



n 



(A5) 



In order to calculate various correlation functions, such as density matrix l|36[) or ground state amplitude l|13|) 
we need to take the limit where several variables t a become equal to each other. This limit is finite, despite the 
apparent singularity in the last integral representation of Z m . We demonstrate this in the simpler case of ground 
state amplitude obtained from l|16|) . Shift the variables t a = t + rj , r\ a — » and rewrite the expression l|A5|l as 



1^] <-'■>■ 



Nm I jtk 



d m xA m {x)- 



fib 

f[ x N e -n*«-«r-ii], 



(A6) 



Due to the presence of totally antisymmetric function A m (x) only the totally antisymmetric part of e tN ^ XaVa survives 
the integration. Using this and the fact that 



,. det a b e lNx ^» N ™(m-l)/2 

n^O A m (irj) lla=o r ( a + 1 ) 

we arrive at the dual representation (|18f) . 

The case of two variables t, t 1 is similar. One starts with the expression l|15|) and shifts the variables 

t a = t + r] a , a = 1, . . . , m/2 

4 = t' + rj b , b = m/2 + 1, . . . , m, 

where r\ a and r\ b go to zero independently. One uses the fact l|Alfl that 

A m (itx, . . .,it m ) = A m /2(ir)a)A m / 2 (iT] b ) [i (t - t')] m /4 
to the leading non- vanishing order in rj a , r\ b and rewrites l|A5|) as 

m/2 r°° A m (xi, . . . ,x m ) e iN ^ XaVa e iN ^ Xb,lb 



(A7) 



Z m (t, t') 



(-iy 



mN I ira 



m/2 



d m x ■ 

[i{t-t')] m2/4 &n{iVa) A m/2 (i Vb ) 

m/2 



a=l 



b=l 



Anti-symmetrising the numerators in each set of variables x a and x b independently and using (|A7|) one gets (1381) . 



APPENDIX B: ANALYTICAL CONTINUATION OF A n 

We use the following integral representation ^(| for the logarithm of Euler's gamma function 

f°° dt /e _2t - e - ' x 

J»r(*) = I 4\ \ - t +(z-i)e- t 

to represent the logarithm of A n as an integral 



t \ l-e~ 



\nA n = V(lnr(a)-lnr(2n + l-a))= / -e^ [ 6 _ J 

a=l Ja 1 L V J 



(Bl) 



(B2) 



where we have summed finite geometric series under the integral. The integral representation defines A n for any value 
of n. In particular for n = 1/2 we get 



In A x 



/2 



dt 



1 — e 2 
1 -e-t 



1 



dt 3e* - 8e*/ 2 + 6 - e"* 



4 7o * (e* - 1)' 



(B3) 
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In order to calculate the last integral we regularise the divergence at t = in the following way: 

hxA 1/a =limC(u) ) C{v)= L - dtt»- l 6e ° e + * (B4) 

4 J ( e * - 1) 

and calculate the integral term by term using the formula 

coo u—l p -p.x J 

r(z/)[C(^-l,/i + 2)-(/i+l)C(^M + 2)], (B5) 



o (e* - if 



where Q(z, q) is the Riemann's Zeta function 

^^E(^j7. C(M)=C0O- (B6) 

The result of integration can be represented as 

C(iz) = 2r(i/) [2(1 - 2"- 2 )C(^ - 1) - (1 - 2^ 1 )C(^) - 1] (B7) 



Taking the limit by the l'Hopitale rule and using the fact that C(Q) = -1/2, CC" 1 ) = -I/ 12 : C'(0) = -ln\/27r we 
arrive at 

InA 1/2 - limC(i/) -3C'(-l) + -^ln2+iln7r, (B8) 
which relates A1/2 to Glaisher's constant A = exp (1/12 — £'(— 1)) as 

Poo = A 2 /2 /\/2 = Tre 1 /^- 1 / 3 ^ 6 . (B9) 

An alternative method of analytical continuation to that described above is to relate the constants A n to Barnes G 
function p|. The definition JSft yields 

4 nLir(a) g(n + l) G(n + 1) G 2 (n + 1) 

" IlLi r (" + a ) G(l) G(2n+1) G(2n + 1)' 1 ' 

where we have used G(l) = 1 and the functional relation G(z + 1) = T(z)G(z). The analytical continuation leads to 
identity A 1/2 = G 2 (3/2) . 
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